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1.0  Summary 

In  this  document,  we  offer  the  final  report  of  the  Defense  Advaneed  Researeh  Projeets  Ageney 
(DARPA)/United  States  Air  Foree  (USAF)  award  #  FA8650-10-1-7045.  An  overview  of  the 
project  is  discussed,  and  the  progress  of  the  research  activity  is  assessed,  as  related  to  the 
program  milestones. 
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2.0  Introduction 


Within  this  project,  a  novel  modeling  approach  is  being  defined  for  the  development  of  a 
transformational  technology  based  on  the  microscopic  thermal  management  of  solid-state 
devices.  By  the  full  inclusion  of  the  phonon  dynamics  within  the  framework  of  charge  transport, 
we  are  implementing  microscopic  models  in  our  existing  particle-based  Cellular  Monte  Carlo 
(CMC)  computer-aided  design  (CAD)  tools  [1]  for  the  design  of  electron  devices  where  the 
generation  and  the  transport  of  thermal  energy  are  functionally  coupled  with  their  electrical 
specifications.  The  crucial  challenges  of  the  present  phase  of  the  project  are  summarized  in  the 
following  three  tasks: 

1 .  Implementation  of  a  rejection  algorithm  for  the  electron-phonon  scattering  table. 

2.  Realization  of  a  solver  for  the  heat  transport  equation. 

3.  Integration  of  phonons  and  electrons  in  a  unified  particle-based  framework. 


Table  1.  Main  Milestones  of  the  Project 


Task 

Milestone 

Deliverable 

Months 
after  award 
(MAA) 

1 

Code  design,  derivation  of 
temperature  -dependent  phonon 
rates 

2MAA 

1 

Implementation  of  rejection  in 
electron-phonon  scattering 

Report  describing  new 
algorithms 

4  MAA 

2 

Implementation  of  a  solver  for  the 
heat  transport  equation. 

Study  of  thermal  conductivity  in 
simple  2D  structure 

5  MAA 

3 

Integration  of  phonons  and  electrons 
in  a  unified  particle-based 
framework. 

Report 

6  MAA 

Accordingly  to  the  statement  of  work,  such  tasks  have  been  organized  in  the  following 
milestones: 

This  is  the  final  report  that  summarizes  the  methodology  developed  to  model  phonon-phonon 
and  phonon-electron  interaction  within  a  perturbative  full-band  approach. 
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3.0  Methods,  Assumptions,  and  Procedures 


3.1  Rejection  Algorithm  for  Scattering 

The  capabilities  of  the  original  CMC  algorithm  [1]  have  been  extended  using  a  rejection 
algorithm.  This  approach  retains  the  simulation  speed  advantages  of  the  CMC  scheme  while 
allowing  the  adaptation  of  the  scattering  rates  to  the  real  time  local  conditions.  Figure  1  shows 
the  rejection  algorithm  flow  chart.  After  the  CMC  triggers  an  event  with  probability  Pcmc  ?  ^ 
scattering  probability  is  computed  from  the  local  conditions.  The  rejection  probability 
Prej^Pioc/PcMc  is  then  used  in  a  stochastic  procedure  to  decide  whether  the  scattering  event 
occurs  or  is  rejected. 

The  scattering  probability  after  the  rejection  is 

Pq^q,  =  Pcmc  *  (^)  =  Pioc  (1) 

where  q  is  the  initial  state,  and  q'  the  final  state. 


Figure  1:  Flow  Chart  of  the  Rejection  Algorithm 
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3.2  Phonon-Phonon  Scattering  Rate 

Perturbation  theory  is  used  to  implement  both  phonon-phonon  and  phono n-defeet  interactions. 
The  anharmonic  decay  and  recombination  scattering  rate  based  on  the  Klemens  model  [2]  is 
given  by  the  following  expression: 


q^q'.q" 


2h  c^iq,  q',  q") 
0)0)  'o) " 


S(a),  0)",  a)’)S(q,  q’,  qf")Fic 


(2) 


where  Fjoc  =  F(n(oj,,  nj"  (,  )  is  a  local  population  dependent  factor,  n(oj,  and  nj"  are  the 
populations  at  q’  and  q",  respectively,  M  is  the  average  mass,  c^(q,  q',  q")  is  the  anharmonic 
coefficient  obtained  from  the  perturbation  Hamiltonian,  and  o),  o)  'and  o)  "are  the  angular 
frequencies  corresponding  to  the  states  q,  q\  and  q",  respectively. 


The  phonon-defect  scattering  probability  based  on  the  work  of  Srivastava  [3]  is  given  by: 


VwF 


'q^q' 


moj'(e*  ■  e^>yS(a)  -  m')(n(oc  +  1) 


(3) 


where  e  is  the  polarization  vector,  is  the  Wigner-Seitz  cell  volume,  and  F  is  an  expression 
that  depends  on  the  phonon-defect  type  to  be  discussed  later. 

Computing  the  local  scattering  probability  requires  estimating  the  local  phonon  populations 
and  n'l'g^  from  the  ensemble.  This  task  requires  defining  a  sampling  volume  of  the  particle 
position  both  in  real  space  (  Vio^)  and  in  momentum  space  (  Hjoc),  then  counting  the  number  of 
simulated  phonons  rj{£Lioc,  Vioc)  iia  the  volume,  and,  finally,  obtaining  the  phonon  population 
using  the  following  approximation: 


^loc  — 


V(.^i  OCf  Vice) 


^state 


OCf  Vice) 


8n^ 

*  - 

^loc^loc 


where  Ngtate  is  the  number  of  states  in  the  volume. 


(4) 


3.3  Phonons  Initialization  and  Boundary  Conditions 

The  number  of  equilibrium  phonons  within  a  sampling  volume  in  a  crystal  of  volume  is: 

^pho(T,Vioc)  = 

p  q 

where  p  is  the  phonon  mode  index,  p,  T)  is  the  Bose-Einstein  distribution,  and  T  is  the 
temperature.  Typically,  Npi^oiT,  Vig^)  is  an  extremely  large  number,  and  even  within  the 
relatively  small  volume  of  a  modem  transistor  there  are  too  many  particles  to  simulate  individual 
phonons  practically.  For  this  reason,  a  weighting  factor  W  is  used,  and  each  simulated  particle 
represents  W  real  particles.  The  factor  W  is  computed  before  starting  the  simulation  using  the 
formula: 
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w  = 


(6) 


^dev 

^sim 


where  N^gy  is  the  total  number  of  real  phonons  in  the  device  and  iV^j^  the  number  of  super 
particles  to  be  simulated. 

The  particle-based  phonon  Monte  Carlo  code  implements  three  boundary  conditions:  (1)  planar 
interfaces,  (2)  rough  interfaces,  and  (3)  contacts.  The  first  boundary  condition  models  ideal 
surfaces,  in  which  the  incident  particle  is  specularly  reflected.  The  second  boundary  condition 
models  interfaces  between  the  extreme  cases  of  perfectly  specular  reflection,  considered  above 
and  perfectly  diffuse  reflection.  The  model  is  based  on  the  work  of  Fuchs-Sondheimer  [4]  [5] 
and  adapted  to  the  full  band  approach.  Contacts  are  modeled  as  ideal  reservoirs  with  infinite 
thermal  capacity  and  constant  temperature.  The  contact  population  is  replenished  injecting 
phonons  in  the  reservoir  using  a  velocity-weighted  probability  distribution,  similar  to  Pardo  [6]. 
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4.0  Results  and  Discussion 


4,1  Validation  of  the  Phonon  Scattering  Rate 

Tests  have  been  designed  to  verify  the  algorithm’s  capability  of  adapting  to  the  real-time  local 
particle  distribution.  The  first  test  verifies  the  capability  of  the  scattering  algorithm  to  reach  the 
equilibrium  steady  state  distribution,  while  a  second  test  verifies  the  capability  of  reproducing 
experimentally  observed  results  over  a  wide  range  of  temperatures. 

Let’s  consider  an  adiabatic  system  with  total  energy  density  E^.  For  any  initial  energy 
distribution,  the  phonons  will  eventually  reach  a  Bose-Einstein  distribution  with  temperature  T^. 
{{Uq  p,  Tc))  that  satisfies  the  following  energy  vs.  temperature  relationship; 

Ec(Tc)  =  (7) 

^pq 

Figure  2  shows  the  energy  distribution  at  the  initial  and  the  final  state  of  the  simulation.  The 
system  is  initialized  with  energy  density  corresponding  to  a  temperature  of  200  Kelvin  (K), 
Ec(200K),  but  the  initial  distribution  is  set  as  that  corresponding  to  a  temperature  of  300K, 
d-ini(.Q>  p)  (^q,p>  300)  (red  line).  The  system  is  allowed  to  evolve  until  it  reaches  steady  state 
(black  line),  at  which  time  the  energy  density  corresponding  to  the  theoretical  distribution 
{riqp,  200/f)(blue  line)  is  obtained. 


energy  [ev] 

Figure  2:  Simulated  Evolution  of  a  System  Initialized  with  a  Non-equilibrium  Distribution 

of Phonons 

The  system  evolves  to  a  final  distribution  coincident  to  the  analytical  reference  distribution. 
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4.2  Phonon-Isotope  Scattering 

An  isotope  is  an  atomic  variant  of  an  element  with  the  nominal  number  of  protons,  but  a  different 
number  of  neutrons  and  hence  a  different  atomic  mass.  For  example,  and  are 

isotopes  of  Carbon  having  12,  13,  and  14  neutrons,  respectively.  The  isotopes  of  an  element 
show  identical  chemical  behavior  (except  in  the  reaction  speed)  as  well  as  both  crystal  and 
electronic  structure,  however,  the  difference  in  mass  disrupts  the  periodicity  of  the  crystal 
resulting  in  a  large  phonon  scattering  rate.  The  overall  effect  of  this  is  seen  in  a  large  reduction 
of  the  thermal  conductivity  tensor  (fc). 

The  thermal  conductivity  of  the  natural  crystal  the  isotopically  pure  crystal  (fcjso)  can 

differ  by  up  to  one  order  of  magnitude  at  low  temperatures  [7].  At  room  temperature  the  effect  is 
more  ambiguous,  and  even  for  a  material  studied  as  extensively  as  silicon  (Si)  the  reported  /Cj^o 
range  is  between  110%  [8]  [9]  and  160%  [10]  [11]  [12]of  the  k^at  value.  For  these  reasons 
isotope  scattering  cannot  be  overlooked  in  the  simulation  of  real  semiconductor  across  different 
ranges  of  temperature.  Isotope  scattering  is  implemented  as  a  phonon-defect  scatter;  using 
Srivastava’s  model  expressed  in  equation  (3)  where  the  parameter  F  is  given  by: 


I 


and  where  /j  is  the  fraction  of  atom  with  mass  Mj,  and  M  the  average  mass. 

Figure  3  shows  the  effect  on  the  thermal  conductivity  of  Si  of  the  scattering  model  implemented 
in  our  particle  based  code.  The  lines  represent  the  measured  [8]  thermal  conductivity  of  natural 
silicon  (92.2%  Si^^,  4.6%  Si^‘^  3.1%  Si^^)  and  an  isotopically  enriched  sample  of  Si^^. 
The  simulations  without  the  isotope  scattering  follow  the  isotopically  pure  thermal  conductivity, 
while  including  the  isotope  scattering  model  allows  for  accurate  reproduction  of  the 
^iuat  conductivity. 
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Figure  3:  Simulated  and  Experimentally  Determined  Effects  of  Phonon-Isotope  Scattering 
on  Low  Temperature  Thermal  Conductivity:  the  Natural  Silicon 

92.2%  Si^^,  4.6%  3.1%  is  compared  to  the  enriched  sample.  The  agreement  is 

excellent 

4,3  Phonon  Dynamics 

The  ability  of  the  Monte  Carlo  approach  to  reproduce  thermal  transients  has  been  tested  on  a 
thermal  resistor  with  a  length  of  2  micrometers.  The  device  has  two  lateral  contacts  set  to 
temperatures  Th  =  310 K  and  Ti  =  290K,  and  the  device  is  initialized  at  the  temperature  T^. 
During  the  simulation,  the  transient  effective  temperature  is  extracted  from  the  local  energy 
density  by  inverting  the  energy  vs.  temperature  relationship  seen  previously  in  equation  (7). 

The  simulated  evolution  of  the  temperature  can  be  estimated  analytically  by  solving  Pick’s  law 
of  diffusion.  An  approximate  analytical  solution  of  Pick’s  law  may  be  obtained  by  means  of  the 
Laplace  transform; 


Tp{x,  t)  —Ti  +  (Tfj  Ti 


(9) 


where  Tp  is  the  estimated  temperature,  t  the  time,  D  the  diffusion  coefficient,  L  the  length  of  the 
device,  and  erfc  is  the  complementary  error  function.  The  complementary  error  function  erfc  is 
related  to  the  error  function  by  erfc(z)  =  1  —  erf(z),  and  the  error  function  itself  is  defined  as; 


(10) 
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Figure  4  shows  both  the  theoretieal  temperature  distribution  (dashed  lines)  and  the  results  of 
Monte  Carlo  simulation  (solid  lines)  being  averaged  over  50  simulations  (eaeh  one  with  different 
random  number  seed).  The  full  band  phonon  Monte  Carlo  is  seen  to  follow  the  estimation  from 
Fiek’s  law  at  eaeh  time. 


Figure  4:  Comparison  between  the  Theoretical  Temperature  Distribution  (dashed  lines), 

and  the  Monte  Carlo  Simulation  (solid  lines) 

The  simulations  follow  the  analytical  solution  during  the  time  transient. 


4,4  Phonon-Electron  Interaction 

Within  the  CMC  framework  the  seattering  probabilities  from  all  initial  states  to  all  the  possible 
final  states  are  pre-computed  for  a  specifie  temperature  and  stored  in  a  lookup  seattering  table. 
Whenever  a  seattering  event  oceurs,  the  final  state  is  ehosen  from  the  lookup  table.  This 
approaeh  would  require  using  a  different  seattering  table  to  model  different  temperatures.  A  new 
structure  for  the  scattering  table  and  an  enhanced  scattering  algorithm  has  been  designed  to 
model  a  wide  range  of  temperatures  using  only  one  table  while  retaining  the  speed  advantage  of 
the  CMC. 


The  new  structure  uses  a  two  level  approach.  Each  element  of  the  first  level  array  contains  the 
information  used  in  a  typical  Ensemble  Monte-Carlo  (EMC),  e.g.,  the  scattering  rate,  scattering 
mechanism,  phonon  mode  involved,  etc.  Each  element  of  the  EMC  array  is  linked  to  a  second 
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array  (CMC  array),  which  contains  all  the  possible  final  states  and  the  scattering  probability  for 
that  mechanism.  The  new  algorithm  requires  two  steps:  first  the  new  scattering  table  is  pre¬ 
computed  according  to  an  expected  maximum  temperature  (Tmax),  second  a  rejection  technique 
is  used  to  adapt  the  scattering  probability  to  the  local  temperature  (Tloc). 

Based  on  this  scheme,  the  scattering  is  modeled  as  a  five  step  process: 

1 .  The  scattering  table  is  pre-computed  for  the  maximum  expected  temperature  Tmax. 

2.  A  scattering  mechanism  is  chosen  from  the  EMC  array 

3.  A  final  state  is  chosen  from  the  related  CA  array 

4.  A  rejection  probability  is  computed  according  to  the  local  runtime  conditions 

5.  A  stochastic  process  establishes  whether  the  scattering  occurs  or  is  rejected. 

Figure  5  shows  the  electron  drift  velocity  versus  electric  field  in  gallium  arsenide  (GaAs) 
computed  with  the  original  CMC  table  structure  (squares)  at  temperature  T=150K,  and  the  new 
structure  with  the  rejection  technique  (solid  lines)  with  Tmax=400K  and  Tloc=150K.  The 
simulations  overlap  at  low  electric  field,  where  the  CMC  grid  is  extremely  fine.  However,  at  high 
electric  fields  where  the  grid  is  coarser,  the  new  structure  produces  a  more  accurate  final  state 
energy. 
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Figure  5:  Validation  of  the  Rejection  Algorithm  hy  Simulation  of  the  GaAs  Electron  Drift 

Velocity  versus  Electric  Field 

The  velocity  computed  using  the  rejection  algorithm  (red  line)  is  identical  to  the  one  obtained 
with  a  fixed  temperature  table  (dots).  Further  improvement  is  obtained  by  the  implementation  of 

the  energy  correction  (dashed  line) 
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Figure  6  shows  the  drift  velocity  versus  electric  field  in  Si  at  four  different  temperatures.  The 
simulations  employing  the  scattering  table  with  the  rejection  technique  (lines)  are  performed 
using  the  same  pre-computed  Tmax=350K  rejection  table  and  varying  Tloc,  while  the  reference 
simulations  (dots)  use  a  separate  table  for  each  temperature. 
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Figure  6:  Electron  Drift  Velocity  versus  Electric  Field  in  Silicon  at  Different  Temperature 

The  lines  are  all  computed  using  the  same  Tmax=350K  scattering  table  and  the  rejection 
technique,  whereas  the  velocities  represented  by  dots  are  computed  using  fixed  temperature 

tables 


4,5  Beyond  the  Heat  Equation:  Energy  Balance  Equation  Solver 

The  original  research  program  proposed  an  efficient  technique  to  reach  thermal  equilibrium 
within  the  device  layout.  This  technique  was  based  on  a  solution  of  the  Heat  Transport  Equation 
(HTE)  (flux -based)  coupled  self-consistently  with  the  particle-based  electron  dynamics.  Once  the 
electro-thermal  steady-state  conditions  were  reached,  the  temperature  map  supplied  by  the  HTE 
solver  would  have  been  replaced  by  a  population  of  phonons  and  the  particle-based  phonon 
dynamics  simulation  engine  would  have  been  started  in  order  to  solve  transients  and  non¬ 
equilibrium  heat  transport. 

However,  an  approach  based  on  solving  the  Energy  Balance  equation  for  phonons,  in  lieu  of  the 
Heat  Transport  Equation,  was  adopted.  The  Energy  Balance  approach  is  better  suited  to  self- 
consistent  coupling  with  the  electron  dynamics  and  allows  a  higher  degree  of  accuracy  by 
supplying  a  separate  solution  for  each  phonon  mode  (or  group  of  modes).  The  approach 
developed  here  is  based  on  the  energy  balance  equation  for  each  mode  p,  directly  obtained  from 
the  phonon  Boltzmann  Transport  Equation  (BTE); 
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e-p 


P-P 


(11) 


dWf, 

5t 


-V  ■  Fjy  + 


dW^ 


St 


dWf, 

St 


where  W^(r,  t)  =  ^Zk  E^(k)  f^(r,  k,  t)  is  the  ensemble  energy  in  the  volume  D.  of  the  reeiproeal 
spaee,  F^(r,  t)  =  -Zk v(k)E^(k)  f^(r,  k,  t)  is  the  energy  flux,  and  the  two  partial  derivatives  of 
in  the  LHS  of  the  equation  represent  the  rate  of  change  of  the  ensemble  energy  density  due  to 
electron-phonon  and  phonon-phonon  interaction,  respectively. 


Under  steady-state  conditions  the  time  derivative  on  the  LHS  must  be  zero,  and  hence  the  heat 
flux  is  given  by  the  sum  of  electron-phonon  and  phonon-phonon  contributions,  i.e.. 


V-(K^(T,r)VT) 


aw„ 


dt 


-f  ■ 


aw„ 


e-p 


dt 


=  -p„ 


p-p/ 


(12) 


where  a  Fourier  Law  approximation  has  been  used  for  F^(r)  =  —k^(T,  r)VT,  and  the  subscript 
p  denotes  either  the  optical  or  acoustic  phonon  mode,  respectively. 


The  CMC  code  includes  efficient  solvers  for  linear  elliptical  partial  differential  equations 
(PDFs),  and  thus  it  is  desired  to  manipulate  the  energy  balance  equation  in  order  to  re-write  it  in 
the  form  of  an  elliptical  PDE.  The  main  issue  with  such  manipulation  is  the  dependency  of  k^ 
from  the  temperature  and  the  position.  We  therefore  assume  that,  with  respect  to  the  position,  the 
thermal  conductivity  can  be  represented  as  a  piece-wise  function  of  the  temperature,  in  other 
words,  k^  c(T)  is  a  function  of  the  temperature  but  is  not  changing  with  the  position  within  each 
cell  C  of  the  finite  differences  grid.  We  therefore  express  this  restricted  position  dependency 
with  the  index  C  rather  than  via  a  full  functional  dependence  on  the  position  vector  r. 


With  this  in  mind,  we  utilize  the  well-known  Kirchhoff  Transformation  to  define  an  “apparent” 
temperature  0  (T )  ^ 

1  r'^ 

en,c(T)  =  To  +r — k^cWdx  (13) 

vml  Ito 

where  Tq  is  a  reference  temperature,  and  k^  c(To)  the  thermal  conductivity  of  the  phonon  mode  p 
at  the  reference  temperature.  This  allows  us  to  rewrite  the  energy  balance  equation  as  follows: 


v^e 


n,c 


P,(r) 

S-c(To) 


(14) 


which  is  a  linear  Poisson  equation  for  the  temperature  0^,c(T).  This  final  equation  can  be 
efficiently  solved  with  the  multi-grid  solvers  readily  available  in  the  code. 


The  temperature  dependent  thermal  conductivity  is  available  for  many  materials  as  a  power  law 
fit  to  experimental  data  [13].  Using  these  relationships,  the  Kirchhoff  Transformation  is 
computed  within  the  code.  This  computation  is  performed  for  two  groups  of  modes,  the  acoustic 
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and  the  optical,  and  it  is  assumed  that  the  temperature  dependence  affects  only  the  acoustic 
modes.  The  temperature  value  along  with  its  corresponding  acoustic  and  optical  Kirchhoff 
temperatures  are  tabulated  for  each  material  at  the  beginning  of  a  simulation  and  used  for  the 
Kirchhoff  and  inverse  Kirchhoff  transformations  throughout  the  simulation. 

To  validate  the  approach,  we  have  solved  the  above  equation  with  the  forcing  function  set  to 
zero,  i.e.,  within  the  diffusive  regime.  In  this  case  we  expect  the  solution  to  simply  be  linear  in 
between  two  boundary  “thermal  contacts”  where  the  temperature  is  specified.  To  verify  our 
approach,  a  thermal  GaAs  resistor  has  been  simulated  where  a  right  contact  is  initially  set  at 
31  OK,  while  a  left  contact  and  all  intermediate  cells  are  set  to  300K.  The  temperature  map  is 
first  transformed  into  a  Kirchhoff  apparent  temperature,  then  the  PDE  is  solved,  and  finally  the 
inverse  transformation  is  applied  to  return  to  the  “standard”  temperature  variable.  The  resulting 
temperature  map  is  linear  as  expected  as  seen  in  Figure  7. 


Figure  7:  Steady-State  Temperature  Distribution  within  a  Device  with  the  Left  Contact  set 

at  300K  and  the  Right  one  at  310K 

The  solution  is  linear  in  the  region  within  the  contacts  as  expected  in  the  absence  of  a  forcing 

function 
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Next,  a  manner  to  calculate  the  forcing  function  for  the  energy  balance  equation  must  be 
implemented.  The  forcing  function  for  a  phonon  mode  p  can  be  separated  into  individual 
contributions  due  to  electron-phonon  (e-p)  scattering  and  to  phonon-phonon  (p-p)  scattering, 
representing  the  rate  of  change  of  the  phonon  energy  density  for  the  particular  mode, 
respectively.  The  forcing  function  is  then  written  as: 


-P,= 


'dW,. 


dt 


dW,. 


+ 


e-p 


dt 


P-P/ 


(15) 


The  electron-phonon  term  in  the  forcing  function  can  be  computed  by  tracking  the  energy 
exchanged  between  phonons  and  electrons  during  scattering,  while  the  phonon-phonon  term  is 
expressed  using  a  relaxation  time  approximation  [14]: 


dWf, 


dt 


=  -C, 


'  op 


op-ac 


op 


^op-ac 


dW^ 


dt 


=  -Cn 


ac-op 


T  —T 

^ac  ^op 


ac-op 


(16) 

(17) 


where  is  the  volumetric  heat  capacity  the  temperature  of  mode  //,  and  Tj_y  the 

energy  relaxation  time  between  mode  i  and  j.  In  the  above  equations,  the  subscripts  ac  and  op 
denote  the  acoustic  and  optical  modes,  respectively. 

4,6  Boundary  Conditions 

It  was  initially  assumed  that  a  simple  Dirichlet  boundary  condition  [15]  [16]  on  a  thermal  contact 
would  be  sufficient  for  our  simulation  purposes.  In  the  Dirichlet  boundary  condition,  the 
temperature  value  is  simply  set  on  the  boundaries  where  specified  and  not  allowed  to  change 
from  this  prescribed  value,  i.e.,  7),  =  c 

However  we  observed  that,  both  in  the  Monte  Carlo  simulations  as  well  as  the  commercial 
simulator  Synopsys  [17],  that  this  condition  seemed  to  cool  the  system  too  much  with  peak 
temperatures  in  a  biased  device  rising  only  a  few  Kelvin  above  room  temperature.  Other 
boundary  conditions  explored  are  those  of  the  inhomogeneous  Neumann  [15]  [16]  and  the  Robin 
[15]  [16]  boundary  conditions.  It  should  be  noted  that  the  Robin  boundary  condition  is  also 
called  a  generalized  Neumann  condition  within  the  Matlab  [18]  software  PDE  package.  In 
addition,  the  Dirichlet  boundary  condition  is  sometimes  called  a  boundary  condition  of  the  first 
kind,  the  Neumann  boundary  condition  a  boundary  condition  of  the  second  kind,  and  the  Robin 
boundary  condition  a  boundary  condition  of  the  third  kind  [16]. 
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The  Neumann  condition  is  that  of  prescribing  the  heat  flux  that  leaves  the  specified  boundary, 
and  is  expressed  as: 


(18) 


while  the  related  homogenous  Neumann  condition,  for  instance  that  of  an  insulating  surface,  is 
expressed  as: 


(19) 


The  problem  with  this  condition  is  that  it  does  not  produce  a  unique  solution,  but  rather  a  family 
of  solutions  differing  by  an  additive  constant  of  integration.  Hence,  to  effectively  use  this 
condition  the  temperature  must  be  prescribed  at  some  location  within  the  simulation  domain  R, 
i.e.,  there  must  be  at  least  one  point  within  the  device  in  question  at  which  we  have  a  Dirichlet 
boundary  condition.  In  addition  the  elliptic  PDE  with  a  Neumann  boundary  condition: 


=  F  inR 


(20) 


dv 

—  —  q  on 

dn  ^ 


dR 


(21) 


is  also  subject  to  the  following  constraint  [15]: 

-  F  d(x,y)  =  jg^g  dS.  (22) 

The  Robin  boundary  condition  is  expressed  as 

-/t£=h(7'-7'.)  (23) 

where  T  is  the  interior  temperature  at  the  boundary  and  is  the  temperature  of  the  outside 
environment  or,  in  our  case,  the  temperature  of  the  heat  sink  of  the  device.  The  heat  transfer 

coefficient,  h,  has  units  and  can  be  thought  of  as  a  thermal  surface  conductance,  or  the 

inverse  of  a  thermal  surface  resistance.  The  role  of  the  heat  transport  coefficient  is  to  scale  the 
outwardly  directed  flux  based  on  how  well  heat  should  be  dissipated  away  from  the  boundary. 

k- 

The  Robin  condition  yields  a  unique  solution  as  long  as  -  >  0  [15]. 

The  Dirichlet  and  Robin  boundary  conditions  have  been  implemented  and  tested  using  a  simple 
case  with  an  analytic  solution.  In  this  test,  a  uniform  piece  of  GaAs  is  used  with  either  the 
Dirichlet  or  Robin  condition  set  on  the  bottom  plane  of  the  material  and  a  homogenous  Neumann 
condition  set  on  all  other  surfaces  (meaning  that  zero  heat  flux  leaves  those  surfaces). 
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The  governing  differential  equation  is  a  Poisson  Equation  of  the  form  (in  1-D) 


d^T  q 

dy^  K 

whieh  has  a  simple  analytieal  solution  of 

7’=  +C^y  +  C2 


(24) 


(25) 


The  eonstant  ,  where  L  is  the  length  of  the  heat-generating  material,  for  both  the 

Diriehlet  and  Robin  boundary  eonditions.  For  the  Diriehlet  eondition,  the  eonstant  C2  is  given 
by  the  temperature  speeified  at  the  boundary,  C2  =  Ti,,  while  for  the  Robin  eondition  it  is  found 
that  C2=Ta+ 

Various  values  for  the  uniform  heat  generation  rate  q  have  been  used  in  the  Diriehlet  ease,  while 
different  values  of  the  transfer  eoeffieient  h  were  used  in  the  Robin  ease.  Simulations  were  then 
earried  out  eonsidering  the  eases  of  a  temperature  independent  thermal  eonduetivity  (as  for  the 
analytie  solution)  and  that  of  a  temperature  dependent  thermal  eonduetivity. 

The  results  of  these  tests  are  shown  in  Figure  8  and  Figure  9,  respectively.  The  numerical  results 
for  the  temperature  independent  case  are  seen  to  overlap  the  expected  analytic  solutions.  Flsing  a 
temperature-dependent  thermal  conductivity  in  the  simulation  results  in  a  rise  in  the  temperature 

distribution  of  approximately  2%  in  the  case  of  a  Robin  condition  with  h  =2e5  0.08% 

with  h  =le7-^,  0.8%  in  the  case  of  a  Diriehlet  condition  with  q  =5.34el6  and  0.0006  % 

with  q  =5.34el3  — .  This  rise  in  temperature  is  expected  due  to  a  decrease  in  thermal 
conductivity  as  the  temperature  increases.  The  difference  between  the  constant  and  temperature 
dependent-thermal  conductivity  temperature  distributions  also  increases  as  the  temperatures 

become  higher,  e.g.,  a  roughly  2%  difference  for  the  Robin  condition  with  h  —2q5  ^  with 

h  =2e5  where  the  peak  temperature  for  the  constant  case  is  350. 75K  but  only  a  0.08%  change 
w 

with  h  =Iq1  —  where  the  peak  temperature  is  301.27K.  Due  to  the  weak  non-linearity  in  the  k 

vs  T  relationships,  the  accuracy  of  the  simulated  result  in  the  constant  k  case,  and  the  behavior  of 
the  changes  in  the  temperature  distributions  of  the  dependent  k  simulations  we  conclude  that  the 
use  of  our  Kirchhoff  transformation  approach  is  justified. 
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Figure  8:  Results  for  the  Dirichlet  Boundary  Condition 

The  constant  thermal  conductivity  simulation  (red)  is  seen  to  overlap  the  analytic  solution  (blue), 
while  the  simulation  with  a  temperature-dependent  thermal  conductivity  reaches  higher  peak 
temperatures  as  expected  due  the  inverse  relationship  between  thermal  conductivity  and 

temperature 
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Figure  9:  Results  for  the  Rohin  Boundary  Condition 

The  constant  thermal  conductivity  simulation  (red)  overlaps  the  analytic  solution  (blue),  while 
the  simulation  using  a  temperature-dependent  thermal  conductivity  reaches  higher  peak 
temperatures  due  to  the  inverse  relationship  between  thermal  conductivity  and  temperature 
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5.0  Conclusions 


This  project  presented  a  novel  modeling  approach  for  the  development  of  a  technology  based  on 
the  microscopic  thermal  management  of  solid-state  devices.  New  techniques  that  model  the 
phonon  dynamics  within  a  full  spectrum  framework  have  been  presented  along  the  comparison 
between  simulated  and  experimental  results.  A  rejection  technique  has  been  implemented  to 
adapt  the  electron  scattering  rates  to  the  local  conditions,  and  the  thermal  conductivity  has  been 
correctly  reproduced  using  phonon-phonon  and  isotope  scattering  within  a  particle  based 
approach. 

In  addition,  an  energy  balance  approach  derived  from  the  phonon  BTE  was  successfully 
manipulated  into  the  form  of  an  elliptic  PDE  through  the  implementation  of  the  Kirchhoff 
Transformation,  and  agreement  with  analytical  solutions  was  demonstrated.  Boundary 
conditions  relevant  to  thermal  simulations  were  also  explored  resulting  in  the  implementation  of 
the  Robin  boundary  condition. 
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LIST  OF  SYMBOLS,  ABBREVIATIONS,  AND  ACRONYMS 


ACRONYM 

DESCRIPTION 

BTE 

Boltzmann  Transport  Equation 

CAD 

computer-aided  design 

CMC 

Cellular  Monte  Carlo 

EMC 

Ensemble  Monte  Carlo 

GaAs 

gallium  arsenide 

HTE 

Heat  Transport  Equation 

DARPA 

Defense  Advanced  Research  Projects  Agency 

K 

kelvin 

MAA 

months  after  award 

PDE 

partial  differential  equation 

Si 

silicon 

Tloc 

local  temperature 

Tmax 

maximum  temperature 

USAF 

United  States  Air  Force 
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